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1.  Introduction 

High  speed  impact  of  an  object  on  solid  media  is  a complicated 
phenomena  which  is  influenced  by  variables  such  as  material  properties, 
impact  velocity,  projectile  shape,  target  support  position  and 
relative  dimensions  of  the  target  and  projectile.  In  this  research 
work  we  have  concentrated  our  effort  in  the  investigation  of  a class 
of  impact  problems  whose  angle  of  impact  is  so  shallow  that  the 
target  impact  disturbance  might  be  considered  as  a surface  wave 
problem. 

The  above  assumption  has  certain  advantage.  The  theory  for  the 
wave  disturbance  resulting  from  the  application  of  a pressure  dis- 
tribution on  the  fluid  surface  is  well  known.  Lamb  [ 1] . Knowing  the 
shape  of  the  projectile  the  resulting  pressure  at  the  contact  area 
can  be  determined.  Under  these  pressure  the  dynamic  behavior  of  the 
projectile  can  be  computed  by  a step  by  step  approach. 

Our  other  goal  is  to  simulate  this  class  of  impact  problems 
using  a general  purpose  finite  element  program.  The  NASTRAN  program 
was  selected  because  this  general  purpose  program  is  most  widely  used 
by  the  Army  ordnance  engineering  community,  and  the  program  has  the 
capability  of  handling  load  matrices  due  to  nonlinearity. 

Before  we  could  start  to  study  the  pressure  developed  during 
impact  we  had  to  be  certain  that  NASTRAN  can  handle  any  nonlinear 
material  problems. 

Most  of  this  first  year  of  research  on  this  project  was  spent 
in  doing  nonlinear  material  problems  using  NASTRAN.  The  theoretical 
background  for  the  pressure  distribution  developed  during  impact  is 
given  in  the  last  chapter.  Unfortunately  no  "numerical  experiments" 
were  performed  to  verify  the  theory  presented. 
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2.  Nonlinear  Material  Problems  using  NASTRAN 


, J NASTRAN  is  a finite  element  computer  program  for 

structural  analysis  that  is  based  on  the  displacement  approach. 

The  program  includes  a limited  capability  for  the  solution  of 
> nonlinear  elastic-plastic  problems  using  Piecewise  Linear  Analysis. 

In  version  16  of  NASTRAN  the  elements  that  do  have  this  capability 
are  ROD,  TUBE,  BAR,  and  PLATE.  In  this  report  a Direct  Iteration 
Method  is  described  that  allows  the  elastic-plastic  analysis  of 
structures  that  can  be  represented  by  combinations  of  elements 
contained  in  the  NASTRAN  element  library,  such  as  beams,  rods,  shear 
and  twist  panels,  triangular  and  quadinlateral  plates,  conical  and 
toroidal  shells,  solids  of  revolution,  scalar  elements,  three 
dimensional  elements,  isoparametric  elements,  and  constraint  elements. 

L. 

In  the  Piecewise  Linear  Analysis  the  load  is  applied 
in  increments  such  that  the  stiffness  properties  can  be  assumed 
to  be  constant  over  each  increment.  The  stiffness  matrix  for  each 
increment  is  dependent  on  the  current  state  of  stress  in  the 
structural  elements.  The  increments  in  displacements  and  stresses 
are  accumulated  to  produce  the  final,  nonlinear  results.  Since  the 
algorithm  assumes  linearity  between  sequential  loads,  the  results 
will  depend  on  the  user's  choice  of  load  increments.  When  the  user 
selects  large  load  increments  and  the  material  properties  are  changing 
rapidly,  the  results  may  be  unacceptably  inaccurate.  If  small  load 
increments  are  used  when  the  structure  is  nearly  linear  the  solution 
will  be  more  accurate  but  relatively  costly.  Based  on  the  same 
approach  Yang  and  Frederick  C 21  expanded  the  NASTRAN  capability  to 
include  axisymmetric  solid  ring  elements. 


The  Direct  Iteration  Method  provides  a very  simple  but 
approximate  solution  of  an  elastic-plastic  problem.  In  brief,  the 
procedure  is  as  follows.  Assume  that  all  material  is  elastic,  compute 
the  structure  stiffness . matrix,  apply  the  full  load  to  the  structure, 
and  solve  for  displacements  and  stresses.  Compute  strains  in  each 
element.  For  those  elements  strained  beyond  their  elastic  limit. 


r 
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approximately  reduce  their  elastic  stiffness,  and  compute  a different 
structure  stiffness  matrix  based  on  the  modified  elastic  properties. 
Again  apply  the  full  load  and  compute  displacements  and  stresses. 

Stop  iteration  when  source  convergence  criteria  is  satisfied.  The 
final  set  of  displacements  and  stresses  is  the  elastic-plastic  solution 
of  the  structure. 

The  principal  advantage  of  the  Direct  Iteration  Method  is 
its  simplicity.  However,  there  are  several  limitations.  The  method 
may  not  converge  if  plastic  region  is  extensive.  Possible  elastic 
unloading  of  plastically  deformed  material  is  not  accounted  for.  The 
flow  rules  of  plasticity  are  ignored  by  the  assumption  that  the  solution 
is  uniquely  defined  by  whatever  final  load  level  is  reached.  It  is 
known  that  the  final  state  of  stress  and  strain  depends  on  how  the 
stresses  and  strains  develop  during  loading  as  well  as  how  much  load 
is  applied.  Another  approximation  is  introduced  by  the  application 
of  deformation  theory,  which  is  exact  only  if  all  stresses  increase 
proportionally,  that  is,  if  in  reality,  at  any  point  and  throughout 
loading,  the  resulting  stress  distribution  does  not  depend  on  the 
order  of  applied  loadings. 

With  the  Direct  Stiffness  Method  the  strain  energy  in 
each  element  is  always  underestimated,  therefore  a higher  stiffness 
is  assumed  for  each  element.  For  most  small  strain  (in  the  order  of 
10~^)  problems  this  approximation  is  not  critical  as  shall  be  demon- 
strated in  the  examples. 


The  elastic-plastic  nonlinearity  of  a structural  element 
is  defined  by  the  element  material  used.  Any  isotropic  material 
may  be  made  nonlinear  by  including  a stress-strain  table  defining 
its  extension  test  characteristics.  This  stress-strain  table  must 
define  a nondecreasing  sequence  of  both  stresses  and  strains.  Because 
the  stiffness  matrix  for  the  first  iteration  step  uses  the  elastic 
material  coefficients,  the  initial  slope  should  correspond  to  the 
defined  Young's  Modulus,  E. 

The  stress-strain  curve  can  be  purely  elastic  or  elastic- 
plastic.  In  the  later  case,  the  stress  at  which  plastic  yielding 
begins  has  to  be  given.  For  this  field  criterion  the  hypothesis 
of  von  Mises  seems  to  be  most  widely  accepted,  i.e,  that  yielding 
starts  when  the  effective  or  equivalent  stress 


(T  = 


J ( $*)+ ($2  ji) '+  ( o~3  fly)  ( i ) 


reaches  the  yield  point  of  the  material.  In  the  above  equation 
(f|  t 0^  and  (Jj  are  the  principal  stresses. 

Once  the  stress-strain  relationship  is  defined  the 
calculation  may  be  described  by  the  following  steps: 

1)  After  applying  the  full  load  to  the  structure  the 

elastic  stresses  and  strains  are  obtained.  Strain  energies  of 

the  individual  elements  are  then  computed.  In  NASTRAN  version  16 

or  later  versions,  strain  energy  can  be  requested  as  part  of  the  output. 

2)  The  equivalent  strains  of  each  element  can  be  calculated 


from  the  strain  energy 


where  V is  the  volume  of  the  element.  To  identify  the  nonlinear  \ 

elements  these  strains  are  compared  with  respect  to  the  stress-strain 
curve . 

3)  Plastic  elements  are  determined  by  calculating  the 
equivalent  stress  of  the  elements.  In  all  but  the  simplest  elements 
the  amount  of  plastic  action  is  likely  to  vary  from  one  point  in  the 
element  to  another.  It  is  therefore  necessary  to  select  a point 
whose  equivalent  stress  will  represent  the  behavior  of  that  element. 

In  simple  elements,  the  centroid  is  a good  sampling  point.  In 
isoparametric  elements  the  average  of  the  Gauss  points  values  may 
be  chosen.  The  general  rule  is  that  sampling  point  should  be  placed 
where  strains  are  likely  to  be  computed  most  accurately. 

4)  For  those  elements  that  are  nonlinear  a different 
modulus  of  elasticity  has  to  be  determined  from  the  stress-strain 
diagram.  In  a typical  element,  the  first  elastic  solution  establishes 
line  OA.  The  corresponding  strain  6 , 


can  be  computed. 

5)  With  the  new  secant  modulus  fcs  an  elastic  analysis 

is  again  carried  out  to  reach  point  C in  the  stress-strain  diagram. 
If  this  element  is  in  the  plastic  range,  the  Poisson's  ration 
is  to  be  equal  to  0.5  (for  solid  elements,  setl^=  0.499).  This 
is  to  take  into  account  of  the  simplified  flow  rules. 

6)  Step  2 is  repeated. 

7)  If  there  is  plastic  action  step  3 is  repeated 

8)  Step  4 is  repeated,  etc. 

9)  The  calculation  is  continued  until  the  stresses  or  strains 
do  not  change  beyond  certain  acceptable  percentage  points. 


3.  Modifications 


In  step  5)  of  each  iteration  a complete  reassembly  and 
reduction  of  the  structure  stiffness  matrix  is  required.  To  avoid 
this  inefficiency,  all  nonlinearity  may  be  imposed  as  additional 
loadings.  If  Tk3  is  the  original  stiffness  matrix  of  an  element, 
the  reduced  element  matrix  is 

[W=f(k)=(i-^)=(te)-[M  <3> 

where  [ ^ *(|-  =*(k) 

Thus  the  equilibrium  equation  at  any  iteration  step  n may  be  written 
as 


where  £ k 3 is  the  original  assembled  structure  stiffness  matrix 
{ D}n  is  the  displacement  at  iteration  step  n 
f Rj.  is  the  total  applied  load 

^he  assemble<^  structure  stiffness  matrix  based  on 


at  step  n-1 


In  structure  stiffness  matrix  £km.Jw-|  only  the  modified  elements 
are  assembled. 

Another  advantage  of  this  approach  is  that  it  offers 
the  possibility  of  under  relaxation.  Unrelaxation  is  accomplished 
in  the  selection  of  previous  nodal  displacement  £D at  step  n 
in  the  above  equation.  The  displacement  can  be  given  as 


{D}„-,  = 0-13)  + /3  {d},,., 


(5) 


where  o<fi  < / . For  p=l  there  is  no  relaxation.  Any  other 

p value  the  load  due  to  modified  secant  modulus  is  reduced  at  each 
iteration  step.  The  purpose  is  to  ensure  that  the  resulting  stresses 
converge  after  each  iteration.  Convergence  is 


Figure  2:  Convergence  and  Divergence 


(6) 


illustrated  in  the  left  hand  figures.  If  there  is  divergence  as 
shown  in  the  right  hand  figure,  no  solution  can  be  obtained. 


There  is  another  way  to  avoid  divergence  of  the  solution. 

This  is  to  apply  a fraction  of  the  load  first.  After  completing 
the  iteration  process,  modulus  of  elasticity  are  obtained  for  the 
next  increased  loading.  This  will  avoid  reading  a point  far  away 
from  the  stress-strain  curve  at  the  first  iteration  step,  and 
yielding  may  occur  in  elements  which  later  in  the  process  must  be 
unloaded.  It  is  assumed  that  with  partial  loading,  the  equivalent 
stresses  in  some  elements  do  exceed  the  yield  point. 

Partial  loading  is  often  necessary  for  problems  where 
many  elements  are  in  the  plastic  range.  However,  satisfactory  results 
can  be  obtained  even  with  rather  coarse  load  steps.  This  point  is 
demonstrated  in  the  examples  given  in  this  report. 


4.  Examples 

Case  I.  Five-Rod  Truss 

This  example  was  selected  because  Yang  and  Frederick  had  worked 
out  this  problem  using  both  their  incremental  approach  and  the  price- 
wise  incremental  option  available  in  NASTRAN  (Rigid-Format  6) . 

The  truss  is  loaded  with  100,000  pounds  as  shown  in  Figure  3. 

The  elastic-plastic  material  properties  for  each  member  are  given  in 
Figure  4. 

The  forces  in  this  truss  are  statically  determinate,  therefore 
the  stresses  do  not  change  with  different  modulus  of  elasticity.  In 
order  to  plot  the  complete  behavior  of  the  members  at  different  bad 
level,  two  loads  (80,000  and  100,000  lb.)  were  applied  to  the  truss. 
Numerical  results  for  direct  iteration  analysis  of  the  truss  are 
listed  in  Table  1. 


Table  Is 

Direct  Iteration  Analysis  of  Five-Rod  Truss 


Load 

MODULUS  OF 

E1=WE4 

ELASTICITY 

Vertical 
Displacement 
of  Node  5 

Stresses 

0l=O2=a3=O40 

5=06 

80,000 

10(10)6 

40 (10) 7 

10(10)6 

2.66(10)“2 

5.657(10)4 

8 (10) 4 

80,000 

4.88(10)® 

21. 3 (10) 7 

10(10)6 

5. 39 (10) “2 

5. 657  (10) 4 

8(10)4 

100,000 

2.82(10)6 

17 . 4 (10) 7 

10(10)6 

1.12  (10) 

7. 071(10)6 

10(10) 

The 

vertical  displacements 

of  joints  2 and  5 are 

plotted  in 

Figure  5 

. The  points 

used  to  determine 

the  different 

modulus  of 

elasticity  are  indicated  in  Figure  6.  Since  the  stresses  do  not 
change  for  the  same  applied  load  there  is  no  need  for  iteration  in 
this  sample.  There  is  no  difference  between  the  results  obtained 
by  direct  iteration  and  the  results  obtained  by  Yang  & Frederick 
using  incremental  approaches.  The  final  iteration  run  on  the 
computer  is  listed  at  the  back  of  this  report. 


(/O) 


Case  II.  Three-Rod  Truss 


The  simple  Three-Rod  Truss  (see  Figure  7)  is  statically 
indeterminate.  This  problem  is  worked  out  by  Martin  and  Carey  (3) 
using  a nine  step  incremental  approach.  The  Ramberg-Osgood  des- 
cription of  the  uniaxial  stress-strain  curve  was  used  to  represent 
the  elastic-plastic  material  behavior: 


t:  = 


£ 

£ 


+ 


3 (£.\n 

1 e UJ 


where  E = initial  elastic  modulus 
°YP  *=  yield  stress 

n = integer  specifying  strain-hardening  characteristics  of 
the  material.  For  this  truss  problem,  E=  11.6  (10) 6 psi,  °Yp=  43000 
psi,  and  n=ll.  The  stress-strain  curve  is  plotted  in  Figure  8. 

Applying  the  load  of  120000  lb.,  six  iterations  were  used  to 
approximate  the  solution.  The  secant  modulus  of  elasticity  of  each 
iteration  step  is  indicated  in  Figure  8.  The  results  are  listed  in 
Table  2.  It  compares  favorabilly  with  solutions  given  by  the 
incremental  method.  The  final  iteration  run  is  listed  at  the  back 
of  this  report. 

Table  2:  Direct  Iteration  Method  of  Three-Rod  Truss 
SECANT  MODULUS  OF 

ELASTICITY  DISPLACEMENT  STRESSES 


Iteration 

£i=e3 

E2 

V1 

ar°3 

°2 

1 

1. 16 (10) 7 

1. 16 (10) 7 

0.606xl0“2 

35147 

70294 

2 

1. 106 (10) 7 

7 . 59 (10) 6 

0.779xl0“2 

43054 

59113 

3 

1.04 (10) 7 

6 . 26 ( 10) 6 

0.818xl0“2 

45832 

55184 

4 

9 . 24 (10) 6 

5.53 (10) 6 

0.995xl0~2 

45956 

55008 

5 

8 . 55 (10) 6 

4 . 98 ( 10) 6 

1.088xl0“2 

46527 

54200 

6 

8.42(10)** 

4 . 63 (10) 6 

1 . 134xl0“2 

47733 

52495 

Result  from 
Incremental 

1.42xl0“2 

47916 

52192 

Approach  1 2] 

</5) 


Case  III.  Thick-Walled  Cylinder 


The  thick-walled  cylinder  under  constant  internal  pressure  is 
another  numerical  example  calculated  by  Yang  and  Frederick  C2] 
using  the  manual  incremental  approach.  With  the  direct  iteration 
method  80  axisymmetric  rectangular  solid  elements  were  used  to 
model  the  cylinder  (Fig.  9).  The  stress-strain  relationship  of 
the  material  is  given  in  Fig.  10.  Starting  with  modulus  of  elasticity 
equaling  30 (10) 6 psi  the  strain  energies  of  the  elements  are 
computed.  Using  the  formula  of  equations  (2)  the  equivalent  strain 
of  each  element  was  determined.  With  this  strain  energy,  secant 
modulus  of  elasticity  was  obtained  from  the  stress-strain  diagram. 
(Fig.  10).  If  the  strain  was  in  the  plastic  range,  Poisson's  ratio 
is  set  to  be  0.49. 

Two  internal  pressure  loads  54,562  psi  and  110,570  psi  were 
applied.  For  the  lower  pressure  only  one  iteration  step  beyond  the 
original  run  was  made.  The  results  are  shown  in  Figures  11,  12 
and  13.  Figures  14,  15  and  16  show  the  results  of  the  higher  pressure 
case  after  three  iteration  runs.  The  theoretical  solutions  were 
taken  from  Yang  and  Frederick.  The  agreement  is  good  enough  for 
most  engineering  applications. 
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5.  Pressure  due  to  Penetration 


For  this  presentation  of  the  theoretical  approach,  a two- 
dimensional  model  is  employed  because  of  its  simplicity.  The  basic 
idea  is  to  determine  the  pressure  distribution  by  matching  the 
disturbance  profile  with  the  profile  of  the  deformed  projectile. 

Both  these  profiles  are  due  to  the  interval  pressure  in  the  contact 
surface.  To  compute  this  pressure  distribution  the  effects  of  a 
group  of  pressure  element  on  both  the  projectile  and  the  target  area 
had  to  be  determined. 

The  different  pressure  elements  are 


1.  Mt  Semi  -Infinite  Treasure  'band . 


2 Fowoa\zL  ~ Infinite  Pressure  T>and- 


3 Complete  Triangle  Pressure  tlement. 


4.  /Ait  Half-TJinngle  'Pressure  Element. 


5.  Jvwart  Half -Triangle 

Tissue  Element. 


T—r.  ~ — ; — 1 1 


From  the  pressure  distribution  diagram  the  solution  to  the  prob- 
lem can  be  obtained  by  summing  the  disturbance  profiles  resulting  from 
all  the  pressure  elements,  and  equating  this  to  the  profile  of  the 
deformed  projectile: 


^ , A &)  A (3)  Ai4)  A (4)  -2,2) 

FK-  +1 P,  ( S,j  - -ij  ) + Pf  ( 5,  - •»,  ) + $ Sx  = ° 

j”12  for  .i  = i to  . 


Here  the  superscript  a means  that  the  variable  has  been  made  dimen- 
sionless using  £ g and  the  pressure  for  that  element.  The  double 
subscript  ij  refers  to  the  influence  at  the  ith  field  point  due  to 
a scource  element  at  station  j.  A single  subscript  i implies  that 
j=l,  that  is,  the  source  element  is  located  at  the  leading  edge.  The 
superscripts  (1)  through  (5)  refer  to  the  five  pressure  forms 
considered.  There  are  a set  of  n equations  and  n unknown  pressure 
points  which  can  be  solved  by  a standard  matrix  inversion  technique. 

The  deformations  w^  due  to  the  five  pressure  forms  have  to  be 
computed.  This  is  done  by  standard  finite  element  analysis. 

The  target  deformations  ..  due  to  the  pressure  forms  are  given  by 


(28) 


Lamb  1 . 


Results  due  to  the  five  pressure  forms  are  listed  below: 


Aft  Semi-Infinite  Band: 


f 3 5/P.  = JC syi ( Ae)  1-  H c- A, )\z  COS ( A.) -I  j 


where  A = xk  and  k = cr/c 
o o o 

in  which  f is  the  density,  g is  the  acceleration  of  gravity,  5 is 
the  target  disturbance,  c is  the  velocity  of  the  travellling  pressure, 
and  x is  measured  forward  form  the  point  of  discontinuity  in  the 
pressure.  Also  H is  the  Heaviside  step  function  and  f ix  one  of  the 
two  auxiliary  functions  for  the  sine-and  cosine-integrals,  which 
are  defined  by 


1 

3 


(A)  = 


= r c-i 
Jo  Art 


(t) 


dt 


Forward  Semi-Infinite  Band: 


^ = P(jS'/v0  = 5giUA/)f(  Ac) 


Complete  Triangular  Element: 


5 = n S/p.  = l S 3-* + ] 

+ 2^11  ff  V (>l^  Uc) ] 

- Lit  fjll-Az)  + H(-A,){2Sin  U,)-  A,  j]| 


(29) 


with 


u;f  ( a)  d\ 

A,  = (x~2a)k> 

Az  = cxr2a)tc 

Aft  Half-Triangular  Element: 

-[  J-  Ji  jj 

K j£  SJK.C.U) A<0  +H(-Ac)[2u.'S(At)-l^jj 


Forward  Half-Triangular  Element: 


W9S/r- 


2^  {tW  ‘l(lA‘’>+  H(-Aj[2J0KAt)->M] 


-[i^i(A,)  + H(-A,)(2  5ii>(Ai)-A,}J 

- L n SV’(A- u)fl  A^fH( - A.)[2Cos( A.) -i|j J 


Besides  the  above  listed  deformation  profiles  of  the  target 
area  due  to  the  five  pressure  elements  the  deformation  w^  of  the 
projectiles  due  to  the  same  pressure  elements  have  to  be  computed 
too.  For  this  purpose,  standard  finite  element  analysis  is  used  with 
the  pressure  elements  as  load  to  calculate  the  deformation  profile. 


(30) 
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